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Abstract. The inspirals of stellar-mass compact objects into supermassive black holes are some of 
the most exciting sources of gravitational waves for LISA. Detection of these sources using fully 
coherent matched filtering is computationally intractable, so alternative approaches are required. In 
(]]], we proposed a detection method based on searching for significant deviation of power density 
from noise in a time-frequency spectrogram of the LISA data. The performance of the algorithm 
was assessed in [2] using Monte-Carlo simulations on several trial waveforms and approximations 
to the noise statistics. We found that typical extreme mass ratio inspirals (EMRIs) could be detected 
at distances of up to 1-3 Gpc, depending on the source parameters. In this paper, we first give an 
overview of our previous work in JU, Q], and discuss the performance of the method in a broad 
sense. We then introduce a decomposition method for LISA data that decodes LISA'S directional 
sensitivity. This decomposition method could be used to improve the detection efficiency, to extract 
the source waveform, and to help solve the source confusion problem. Our approach to constraining 
EMRI parameters using the output from the time-frequency method will be outlined. 
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1. BACKGROUND 

Astronomical observations indicate that many galaxies host a supermassive black hole 
(SMBH) in their center. The inspirals of stellar-mass compact objects into such SMBHs 
with mass M ~ few x lO 5 A/ -lO 7 A/0 constitute one of the most important gravitational 
wave (GW) sources for the planned space-based GW observatory LISA. Typical EMRI 
events will come from objects on elliptical orbits around spinning black holes and this 
leaves a significant imprint on the emitted GWs. The parameter space of possible GW 
waveforms is large since the masses of the two components, the amplitude and relative 
direction of the SMBH spin, the orbital pericenter and eccentricity and the object's initial 
phase relative to pericenter all need to be specified. 

EMRI waveforms are complicated. There are three characteristic frequency compo- 
nents arising from the orbital motion, the nodal precession due to spin-orbit coupling 
and general relativistic perihelion precession. The GW power is therefore spread over 
a wide frequency range at harmonics of these three major frequencies. Over an ob- 
servation lasting LISA's nominal lifetime of three years, there is also significant time 
evolution of these frequency components and variation in their relative amplitudes (H]. 
Moreover, due to LISA's motion around the sun, the detector's response to GWs under- 
goes a time-varying amplitude modulation and frequency Doppler shift that depends on 



the source direction. The signal-to-noise ratio (SNR) of a "typical" EMRI event (e.g., a 
10+ 1 6 M & system at 1 Gpc) is large if we can perform an optimal search using coher- 
ent matched filtering. However, the signal power in each Fourier frequency bin is very 
small, given the signal's large frequency spread over N~Tf~ 10 5 bins. 

Study has shown that detection of EMRIs using the optimal fully coherent matched 
filtering is computationally impossible [4], so alternative search techniques are required. 
A semi-coherent matched filtering algorithm has been proposal as an alternative [4]. This 
involves matched filtering on short segments of the data, that are as long as is allowed 
using reasonable computational resources. The power in these segments is then added 
incoherently over the entire data stretch. Preliminary results for this search i4\ suggest 
that the LISA EMRI detection rate will most likely be dominated by inspirals of ~ 10 
M Q BHs onto ~ 1O 6 M SMBHs, and could be as high as ~ 1000 in 3-5 years within 
-3.5 Gpc. 

In this paper, we first in Section [2] give an overview of the time-frequency method 
we proposed and studied previously JTL @L We then discuss the performance of the 
method compared to the (semi) coherent search in Section[3l In SectionS we discuss the 
information about an EMRI event that we can extract using the time-frequency method. 
In particular, we introduce a simple decomposition method to decode the directional 
sensitivity. We also outline the information about the dominant frequency components 
and their evolution, and the evolution of GW power, that can be derived from the data 
once a detection is made. We outline possible applications in Section [51 



2. A TIME-FREQUENCY METHOD TO DETECT EMRIS 

In [Q], 0], we proposed an alternative method for EMRI detection based on the incoherent 
summation of power in a time-frequency plane. The t-f power spectrum is produced 
by dividing the data into segments of equal duration and carrying out a Fast Fourier 
Transform (FFT) (or alternative spectral decomposition technique) on each. In the low- 
frequency regime, LISA can be regarded as a network of two Michelson interferometers 
(denoted I and II here), rotated at 45 degrees relative to one another qQ. The power 
spectrum of the detector is defined for each time segment i and frequency bin k as, 

P(j, k )=2^- + 2 1 -^-, (1) 

°Ik °IIk 

where d{' J denotes the Fourier amplitude of the 7-th segment of data from the 7-th 
detector, of k is the expected variance of the noise, n ! k , in the 7-th detector at frequency 
bin k, assuming the noise is stationary and Gaussian. Our noise model includes the 
unresolvable background from white dwarf - white dwarf binaries in the usual way 
In CO, the powers of the two data streams have been added directly. This is optimal based 
on the maximum likelihood ratio if we have no information about the detector's response 
and simply assume that the signals included in these two data streams are statistically 
independent from one other. 



The strategy is then to calculate the power "density", p (z, k) , by computing the average 
power within a rectangular box centered at each point (z, k), 



n/2 1/2 

P(i,k)= £ £ P(i + a,k + b)/m, 



(2) 



a=-n/2 b=-l/2 



where n, I are the lengths of the box in the time and frequency dimension respectively 
and m = n x / is the number of data points contained in the box. To search for a possible 
signal, we vary the size of the box, and for each choice of n and Z, we search for 
any points at which p (i, k) exceeds a threshold determined by a specified false alarm 
probability set equal for all boxsizes (i.e., we assume an equal chance for any particular 
clustering of the signal in time and frequency). A detection occurs when one such event 
happens in any one box size, (see |0] for details). This approach is very similar to the 
'excess power' technique used in LIGO data analysis [g], which was designed to search 
for significant clustering of excess power in the data caused by a source of unknown 
waveform, but in a given window of time and frequency. Our method gives a simple 
estimate of the power density at each point of the time-frequency plane and therefore 
helps to trace the signal power along the source trajectory. 



3. PERFORMANCE OF THE TIME-FREQUENCY METHOD 



The performance of the time-frequency method can be estimated as follows. When 
the time-frequency window is "wrapped tightly around" most of the signal power, the 
accumulated signal power increases in proportion to m, the size of the window, but the 
noise power is a % 2 distribution with Am degree of freedom, the mean and standard 
deviation of which are 4m and \f%m respectively. The SNR therefore scales with y/m 
and can be written as 



where p 2 is the optimal SNR 2 if we performed matched filtering on the signal within 
the time-frequency window in consideration. 

The required SNR for a detection via this time-frequency method compared to the 
fully coherent method can be understood as follows. In the absence of a signal, the 
output of a matched filter is a Gaussian with zero mean and variance p m . The presence 
of a signal increases the mean to p, 2 and therefore the SNR is p m . If m is large, the 
output of the time-frequency search in equation [2] in the absence of signal will also 
be approximately Gaussian, with mean 4 and variance 8/m. The presence of a signal 
enhances the mean by p, 2 /m. Thus, for a given FAR the optimal SNR required for 
detection by the time-frequency search, pjf, and the optimal SNR required for detection 
using fully coherent matched filtering, p c , are related by 
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(3) 



Ptf (8m) 
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for a given FAP. 



(4) 



The m values for a typical source are around 200-1000. We note that this is under the 
approximation of a large window size and assuming we have included roughly all of 
the signal power. We have also not taken account of the number of trials needed to 
find the source in either case. In the limit where a signal is uniformly distributed over 
the whole time-frequency plane, the performance of the time-frequency method is at 
its worst compared to the coherent method. However, in another limit, not reflected in 
the preceding equation, for a monochromatic source, the time-frequency method can 
perform almost as well as the coherent method (provided we somehow know to use a 
box with m = 1, and that the coherent method still needs to use the power at the end 
of the data stream for detection). In summary, the time-frequency method works better 
when the signal power is well localized. The actual performance depends on how well 
we can find a window that encloses most of the signal power to maximize the SNR. 

The detection efficiency of the time-frequency method has been studied in detail in our 
previous work [2], where only rectangular t-f windows were used for simplicity. We used 
in our simulation a wide range of possible EMRI signals from the 'numerical kludge' 
approximation [3, [a @]. We find that this algorithm is able to detect many different 
EMRI events out to distances of 1-3 Gpc, depending on the source parameters. In an 
untargeted search, a typical source can be detected at 2 Gpc with a detection rate of 
60% at an overall false alarm probability of 1%. Lower eccentricity sources, which have 
less frequency spreading, can be detected as far away as 3 Gpc with a detection rate of 
50% at the same overall FAR The reach of the search can be extended by increasing the 
allowed FAP or by using a targeted search. By comparison, the semi-coherent matched 
filtering algorithm [0] can reach ~ 4.5 Gpc for an overall FAP of 1%, but at a presently 
undetermined detection rate (perhaps ~ 50%). Note that the performance of the time- 
frequency method using a rectangular window is limited by the evolution time scale of 
the EMRI event. 

Broadly speaking, this time frequency search with a rectangular window has better 
than half the reach of the semi-coherent search [4], but at a tiny fraction of the com- 
putational cost. The method can be optimized from using rectangular window by e.g., 
following the trajectory. Many image processing methods can be used to improve this 
shortcoming which we will not discuss in this paper (one possibility is the Hierarchical 
Algorithm for Clusters and Ridges, which will be described elsewhere ifiolp. Given the 
simplicity of the technique, we argue that time-frequency methods could be a valuable 
first step for detecting the loudest EMRI events in the LISA data stream. 

4. EXTRACTING INFORMATION FROM THE 
TIME-FREQUENCY SPECTROGRAM 

There are several issues we must consider when we use the time-frequency method to 
analyze realistic LISA data. One major challenge is to use the results from the time- 
frequency method to extract informations about the source which can then be used to 
reduce the parameter space for a follow up analysis. The results can help reduce the 
computational cost for a subsequent coherent or semi-coherent search. 

In the following subsections, we will first introduce a simple decomposition method 
to decode the directional information for EMRIs. We will then discuss briefly how the 



outputs can be used to optimize the detection statistic, extract wave information and, to 
some extent, resolve overlapping sources from distinct directions. We will emphasize 
EMRI parameter extraction using the t-f method. 



4.1. Directional Information 

LISA can be considered as a network of two detectors in the low-frequency limit 
and three at high frequencies. Therefore, we make use of a network analysis carried out 
for ground-based GW detectors flUl . This is summarized as follows. The response of a 
network of detectors to a source in a given direction is linear. The response in the k-th 
frequency bin and 7-th segment of time can be written in the frequency domain, after 
applying appropriate time delays, as 

dl=A J k h J k + n k , (5) 
For the low frequency 2-detector model of LISA at described in we define 

d{ k /oik\ A j_( f( + /o lk fi x /a lk \ y_( h\ k \ f n\/c lk 

dk ~ \ d{Ja 2k )' k ~\ f{+/c 2k fl*/** J ' * " V K J ' V 4/o 2k 

(6) 

where f- +x are the antenna beam patterns of the i-th detector to the h J +k and h J xk 
polarizations of a GW from a given direction, at time segment j. For more than two 
detectors, the same formalism can be generalized by adding rows in the data and 
noise vectors and in the response matrix. Using a singular value decomposition, the 
response matrix A J can be written as A 7 = Msv r lll2ll where u and v are unitary matrices 
(uu T = /, vv T = I), and s is a diagonal matrix with diagonal values S{ > s 2 called singular 
values. The signal power from equation Q] can then be rewritten as a summation of two 
terms, 

P(j,k) = \(u T dl) l \ 2 + \(u T dl) 2 \ 2 (7) 

with, 

(u T d J k ) 1 = s[h\ + (u T n k ) 1 , (u T d J k ) 2 = s J 2 h 2 + (u T n k ) 2l (8) 

where h\ 2 are the two components of the h' = v T h k . We note that each component of 

h! can be measured statistically independently. Each term in P(j, k) includes a signal 
and a noise with the signal proportional to the square of the singular values s\ 2 , and a 
noise term of unity variance. The singular values s\,s 2 for a given source direction are 
plotted in Figure [Qfor the LISA configuration described in [5]. It is clear that for most of 
LISA's orbit (regions of red-color), the two directional sensitivities are comparable and 
therefore the power summation in equation © is optimal. However, for a significant 
fraction of time (region of blue color), one sensitivity can be much less than the other 
and in some parts of the time-frequency plane, the sensitivity can be nearly zero. 

The significance of this is that the summation method in equation \T\ can be further 
optimized by summing up only terms with large singular values to optimize the SNRs. 
In other words, we use the data from one or both of the time-frequency planes shown in 



Figured] only when they have good sensitivity to the sky direction under consideration. 
It is clear that the improvement in detection efficiency for this simple 2-detector approx- 
imation is limited. However, the situation will be different for configurations accounting 
for the rotation of the LISA constellation (which encodes more directional information), 
and at higher frequencies when LISA becomes effectively a 3-detector network. For a 
3-detector network, at least one null-stream (a term with zero singular values) that is a 
linear combination of the three data streams can be constructed. The detection statistics 
should improve more significantly in this case. The null-stream can also be used to lo- 
calize the source and as a consistency check to run in parallel with the detection method 

fil. 

A direction specific search is required for such a decomposition to work. Based on 
Figure Q] and previous studies on the lower limit of LISA'S angular resolution Q3[| and 
work by Pai et al.yjfl, the angular resolution of LISA will be quite rough for a time- 
frequency method. A search grid to cover all sky directions can be as coarse as tens 
of degrees and the search over multiple directions is therefore not too computationally 
intensive. 



4.2. Constraints on Parameter Space 

Three types of information can be obtained from the time-frequency method once a 
detection is made — (1) f n {t), the time evolution of the dominant frequency components 
in the time-frequency plane, where n indicates the number of harmonic this frequency 
is of its fundamental one. An example of this can be seen in the bright "typical" case 
described in our previous work QjJ, |2|] . (2) f„(t), the time derivatives of the dominant 
frequencies. This can be calculated, e.g., by finite differencing. (3) < \h(t,f„)\ 2 >, the 
signal power in the dominant frequencies, where the brackets <> indicate the average 
over the time-frequency windows that yield SNRs above the threshold and along the 
trajectory. The signal power \h(t,f n )\ can be extracted by inversion of equation ([8]). 
The measurement uncertainty can be calculated based on the singular values. Small 
singular values should be discarded using the standard treatment for an ill-conditioned 



matrix 1112(1 . The continuous signal power, < \h(t,f n )\ >, can then be estimated from 
its discrete representation. 

A full analytical solution to the dynamical evolution of these three quantities is not yet 
available for EMRIs. We will therefore make use of the post-Newtonian(PN) treatment 
described in [3] to illustrate how the time-frequency method might be used to put 
constraints on the system parameters. First we must understand the dependence of the 
dominant frequency components, f n (t), on the system parameters. For a circular orbit 
around a non- spinning black hole, the dominant GW power output is at twice the orbital 
frequency. However, when the central black hole has spin and the orbit is eccentric, the 
GW output is significantly colored by precession effects. The majority of the GW energy 



is emitted close to periapse, when the object is on a nearly circular "whirl" orbit H141I15I1 . 
The dominant GW emission will therefore be at the effective "whirl" frequency and its 
harmonics. We are currently developing expressions for the dependence of the whirl 



frequency on system parameters. In the PN formalism [3], we can write 



f n (t)=nv(t) + f(t)/7t = F l (v(t),e(t),n,M,S,X), (9) 

where v(t) is the radial orbital frequency at time t, and f(t) is the rate of pericenter 
precession (given by equation (29) of |30). The value of n that determines the dom- 
inant frequency component (i.e., the "whirl" frequency) depends non-trivially on the 
eccentricity and spin. The dominant frequencies, f n (t), are therefore functions of v(t), 
eccentricity e(t), reduced mass /i, total mass M, magnitude of the SMBH spin S, and 
X, the inclination of the orbit with respect to the spin direction of the SMBH (assumed 
fixed in H). 

The time derivatives, f n (t), can also be related to those of the system parameters 

fn(t)=F 2 (v(t),e{t),n,M,S,X), (10) 

where, 

dF\ dF\ i d F\ 
F 2 = v^±+e^- + X—±. (11) 
dv de dX 

The time derivatives v, X and e are results of gravitational radiation reaction. PN 
expressions for them are given in equations (28)-(30) of [3] where X is assumed to 
be fixed. 

A relation between a GW frequency, its derivative and the signal power \h(t,f n )\ 2 at 
that frequency is 

< \h(tj n )\ 2 >= ^^E n /(f n f 2 ) =F 3 (v(t),e(t),^M,S,X,D). (12) 

where D is the distance to the source, Ft, = G/ (n 2 c 2 D 2 )E n / (F1F2). Under the quadrupole 
approximation, the energy power radiated at the n-th harmonics can be written as iflol 

32 C^fi 

E n = -^M 2 ^ 4/3 (2^v(0) 10/3 s(H,e(0). (13) 

As an illustration, we show in Fig. [2] the time-frequency power density for a "typical" 
EMRI source of 10+ 10 6 M© at d = 0.5 Gpc (same as Fig. 1, left panel of (^\). Trajecto- 
ries of frequency evolution vs time and harmonic structures are apparent. The weighted 
average of the quantity f 2 \h k \ 2 for the same source are then calculated for each point 
on the t-f plane using Eq. [81 assuming that we have obtained directional information. In 
Fig- El we plot f 2 \h\^ vs frequencies at four different times. The signal power distribu- 
tions among frequency harmonics and their evolution are clearly observable. Data points 
of signal frequencies, their derivatives (f n , f n ) and the fundamental frequency v, can be 
easily measured from the figure. Once the harmonic structures are identified, eccentric- 
ities can be estimated from the relative power in different harmonics (from function 
g(n,e) in Eq.QSfr and fiM 2 / 3 /D can be estimated. A quantitative study for this approach 
is under way. 

In summary equations ©, (flOl) and (PT21) provide an illustration of how constraints on 
the parameter space and information on e(t),v(t) etc. can be derived from a t-f map. 
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FIGURE 1. Singular values s\ (upper panel) and S2 (lower panel) of the response matrix A (Eq.[8]) in 
the time-frequency plane. The EMRJ source is placed at longitude (j) = 60", and latitude 9 = 51", Ecliptic 
orbits. LISA'S antenna beam pattern is taken from [5]. Higher singular values (redder color) indicate higher 
directional sensitivities of LISA. 




FIGURE 2. The t-f power density for a typical EMRI source at d = 0.5 Gpc (see text). Harmonic 
structures of signal frequencies and yearly amplitude modulation caused by LISA'S directional sensitivity 
are evident (cf Fig.[T]). Optimistically, we expect <^ 3 such events in three years. 



If we assume that we have detected Nh harmonics and that along each trajectory we 
have N t independent data points above the detection threshold, then in principle 3N t Nh 
data points can be used to constrain the 2./V + 5 unknowns v(Y),e(?),/l,M, 5, A, and D 
If the SNR is high enough, these parameters can then be calculated using, e.g., the least 
squares method. 
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FIGURE 3. The normalized energy power f 2 \hf\ 2 vs frequency at four different times for the same 
EMRI source of Fig. [2] at d — 0.5 Gpc. The 3rd -6th harmonics of the radial frequency can be identified 
at t = 1.71 yr (solid curve) and at t — 1.84 yr (dashdot line). The four peaks at later times of t = 2.1 T yr 
(dotted line) and of t = 2.81 yr (dashed line) are very close to the 2-5th harmonics of the radial frequency. 



5. DISCUSSION 

We have reviewed a time-frequency method that can be used to detect bright EMRI 
sources. We have discussed the performance of this method — previous studies have 
shown that this method can be effective in detection of EMRIs at distances of up to 
1-3 Gpc depending on the system parameters. We have proposed a recipe based on 
singular value decomposition to decode the directional information about EMRI sources. 
We have presented a method to constrain the parameter space using the PN evolution 
equations given in yD. The parameter space can be constrained using equations ©, 
(flOl) . and (fT2i For a detection of Nf, harmonics with N t data points on each trajectory 
(assuming all trajectories are generated by the same EMRI), there are in principle 3iVyVj 
data points that can be used to constrain the 2N t + 5 unknowns e(t),v(t),H,M,S,D and 
A for the PN expressions in [3]. Therefore, in this simplified model, it is possible, 
in principle, that the total number of observed data points from the time-frequency 
method can be larger than the number of unknowns and the equations can be solved. 
A quantitative study is underway and results will be presented in a follow-up paper. 

Another unsolved issue is the problem of confusion caused by the WD-WD binaries 
in the LISA data. Three different types of information obtained from the time-frequency 
method can be used to distinguish them from EMRIs. One is the frequency binsize 
which gives the maximum SNR, the other is the shape of trajectories and the third 
is the directional information. The maximum SNR of a signal can be obtained only 
in the correct direction and with box sizes smaller than the signal spread. WD-WD 
binary signals are expected to be single-bin (since the Doppler shift is negligible for 
our frequency bin sizes) and the trajectories should exhibit very little time evolution. 

At the end, the time-frequency method provides directly measurements of EMRI in 



dominant frequencies and their time derivatives. These information can be used to map 
out the space-time geometry of the SMBHs lll7L ll8h. This will be the subject of different 
studies (e.g., lH3l). 
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